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A new algorithm is presented, which allows to calculate numerically the partition function Z q of the 
d-dimensional q-state Potts models for arbitrary real values q > at any given temperature T with 
high precision. The basic idea is to measure the distribution of the number of connected components 
in the corresponding Fortuin-Kasteleyn representation and to compare with the distribution of 
the case q = 1 (graph percolation), where the exact result Z\ = 1 is known. As application, 
d — 2 and d — 3-dimensional ferromagnetic Potts models are studied, and the critical values q c , 
where the transition changes from second to first order, are determined. Large systems of sizes 
N = 1000 2 respectively N = 100 3 are treated. The critical value q c (d = 2) = 4 is confirmed and 
q c (d = 3) = 2.35(5) is found. 
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The partition function is a quantity of fundamental 
importance, because it describes completely the behavior 
of any statistical physics model. Unfortunately, in finite- 
dimensions, only few models are analytically tractable 
0|. Hence, Monte Carlo (MC) simulations 0,13 are usu- 
ally applied. The standard approach to obtain the par- 
tition function is to measure the free energy by thermo- 
dynamic integration of the specific heat, i.e. the fluctua- 
tions of the energy. Since this approach is based on mea- 
suring fluctuations, it is not very efficient, hence limited 
to small sizes. One can speedup simulations for certain 



types of systems by applying cluster algorithms 
multi- histogram methods 7] , multicanonical simulations 
0, or transition-matrix Monte Carlo ^| , but the gen- 
eral problem of the strong fluctuations remains. To over- 
come this problem, recently Wang and Landau intro- 
duced 01 a simple yet very efficient method to obtain 
the partition function. The key idea is to measure the 
density of states by performing a biased random walk 
in energy space via spin flips, ft works well for unfrus- 
trated systems, e.g. the standard g-state Potts model 
which has become a standard testing ground for 
Monte Carlo algorithms. The Potts model is of profound 
interest, because, for dimensions d larger than one, it 
exhibits order-disorder phase transitions |13l |. which are 
of second order for q smaller than a critical value q c (d), 
while they are of first order for q > q c (d). ft is analyt- 
ically proven 0] that q c {2) = 4, but e.g. for d = 3, 
the exact value of q r is not known. From various ana- 
lytical work llollIH.Il7| and simulations of moderate-size 
systems [Tsl lia I20j7~2 < g c (3) < 3 seems likely. Un- 
fortunately, since the Wang-Landau method is based on 
spin flips, it works only for integer values of q, hence the 
partition function for 2 < q < 3 cannot be obtained for 
large systems in this way. 

In this letter, an algorithm is presented, which allows 
to calculate numerically the partition function Z q of the 



(i-dimensional g-state Potts models for arbitrary real val- 
ues q > at any given temperature T with high pre- 
cision for large systems. The basic idea is to measure 
the distribution of the number of connected components 
in the corresponding Fortuin-Kasteleyn (FK) represen- 
tation [2l| and to compare with the distribution of the 
case q = I (graph percolation), where the exact result 
Z\ = 1 is known. Large system like N — I000 2 resp. 
N = 100 3 can be treated here, because for the MC simu- 
lation, the cluster algorithm of Chayes and Machta is 
applied. Using this combined approach the critical value 
q c (d = 2) = 4 is confirmed and q c (d = 3) = 2.35(5) is 
found. The outline of the paper is as follows. Next, the 
model is defined. Then the algorithm for calculating the 
partition function is presented. In the main part, the re- 
sults for the d = 2 and d = 3 Potts models are shown. 
Finally a summary is given. 

The g-state Potts model 0] for integer values of q 
consists of N spins Ui € {I, . . . , q} living on the sites of 
an arbitrary graph or lattice G, with the Hamiltonian 
H = — j) &<?i-<?j ) where the sum runs over the edges 
of G, and 8 is the Kronecker delta. For G, here d- 
dimcnsional hypercubic lattices having periodic bound- 
ary conditions with nearest-neighbor interactions are 
considered. The partition function Z q = ^2{ a .j e~ H ^ T , 
T being the temperature, can be written in the FK rep- 
resentation |21| as 
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= P Nb{G ">(l ~ p )N b (G)-N b (G') q Nc(G') ^ ^ 
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where the sum runs over all subgraphs of G having the 
same set of sites and any subset of edges, W q {G') is the 
weight of graph G',p=l- e' 1 / 7 ', N b (G) resp. N b (G') 
are the number of edges in G resp. G' and N C (G') is 
the number of connected components in G' . The FK 
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representation allows for an extension of the model to 
arbitrary real values of q > 0. 

For 5 = 1, Wi{G') is the probability of the subgraph 
G", if the graph is generated randomly by making every 
edge a member of the subgraph with probability p. This 
allows for a very efficient generation of graphs distributed 
according to Wi, i.e. by importance sampling. Since one 
generates each time with probability one some random 
subgraph in this way, one has trivially Z\ = l. 

This allows for a calculation of the partition function 
Z q for any q > in the following way. Let P q (c) the 
probability to have c connected components in a sub- 
graph generated according to weight W q . Then we have 
by definition 
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Hence, we get 



Z q =q 



P q (c) 
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This means, by measuring the probability distributions 
of the number of connected components for random sub- 
graphs (q = 1) and for the target value q, we can obtain 
Z q . Note that Eq. holds for all values of c simulta- 
neously [23. Therefore, by comparison of the full distri- 
butions, one has a mean to determine Z q with very high 
precision . 

Eq. might be useful for analytical calculations, but 
for most interesting graphs G, the distributions cannot 
be obtained in this way. Hence, one uses numerical sim- 
ulations to obtain the distributions Pi(c) and P q {c). In 
practice, one can generate random subgraphs according 
to W\ with importance sampling, as explained above, 
and according to W q using the very efficient cluster algo- 
rithm of Chayes and Machta [5j. This algorithm allows 
for simulation for arbitrary values of q, similar to other 
approaches [H E3 . 

Nevertheless, for large values of q and finite statistics, 
P\{c) and P q (c) will not overlap, because deviations from 
the typical value are exponentially suppressed. In this 
case one has to study intermediate values q\ , . . . , qt £ 
[l,q], calculate each time PqAc) and Z qi . This allows to 
extend P\ (c) stepwise jljj for larger values of c, until 
P\{c) and P q (c) have sufficient overlap. In principle, it is 
a bit ugly, that one has to perform simulations at several 
values of , but on the other hand, one gets the partition 



function for all considered values, which will be useful in 
the following. Note that for the Wang-Landau algorithm, 
also one long run is sufficient only in theory, in practice, 
if the system size is larger than tiny, one has to divide the 
energy range into intervals, perform independent runs for 
each interval, and match the results of the different runs 
as well. Anyway, this is no problem for cither method, 
because it can be done automatically by a program, no 
matter how many intervals have to be matched. The real 
advantage of the present approach is that it works for all 
values of q > 0, since it does not rely on flips of spins. 



-1.9 



T 



T 



T 



1000x1000 Ising model 



exact 

numerically 




FIG. 1: Free energy F per spin as a function of the temper- 
ature T of the two-dimensional Ising model (g = 2) obtained 
by the algorithm and by an exact calculation for system size 
L = 1000. The inset shows the relative error e(T). For each 
independent value of temperature, only a total of 5.5 x 10° 
MC sweeps were performed. 

To test the new approach, it is now applied to the 
two-dimensional Ising model (q — 2), where exact re- 
sults are available for finite system sizes j^. In Fig. 
n the Gibbs free energy per spin F/N = —jjlnZ q is 
shown in the Ising representation (i.e. for the Hamilto- 
nian H = — ■^[28 ai rJj — 1]). The data of the simula- 
tion and the analytical result are given for a large system 
size N = 1000 x 1000. Thus, k = 110 different values q t 
are necessary for measuring P(c) over the desired range. 
Equilibration of the cluster MC simulation is determined 
by monitoring the number of connected components and 
the number of edges when starting with a full resp. empty 
subgraph. Equilibration is assumed, when the values for 
the different starting conditions agree within the range 
of fluctuations. Due to the global update nature of the 
Chayes-Machta algorithm, this is the case for typically 
few Monte Carlo sweeps. Hence, for each value of q t , 
5 x 10 3 steps where sufficient, to obtain a high accuracy, 
as shown in the inset of Fig. ^ 
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FIG. 2: Free energy per spin of the two-dimensional Potts 
model (3 < q < 5) for system size L = 1000. The inset shows 
the distribution of the number of edges at the temperatures 
T(q = 4.00) = 0.910289 resp. T(q = 4.05) = 0.906865. 

Since the aim is to determine q c , the new approach 
is further tested by applying it to the two-dimensional 
Potts model, where q c = 4 is known. In Fig. |3 the 
Gibbs free energy per spin is shown for values in the 
range 3 < q < 5. Due to the large system size, in total 
k = 261 different values for q 6 [1,5] are necessary. For 
q > 4 a kink at the transition temperatures is visible, 
as expected. One could in principle take derivatives of 
the free energy to obtain e.g. average energy and specific 
heat, e.g. to calculate critical exponents in the case of 
second order transitions. A better way is to calculate 
these derivatives analytically from Eq. Q), which allows 
to express the mean energy resp. the specific heat by 
the average number of edges resp. the fluctuations of the 
number of edges 0, 0] , which are available directly from 
the simulation. Since this is a standard approach, it is 
not further pursued here. 

One can determine q c more precisely by considering 
the distribution of the number of edges fl8t H^ . Ejf , see 
inset of Fig. [3 The distributions are obtained by per- 
forming several long simulations for T G [0.906,0.907] 
resp. T £ [0.9100, 0.9105] , exhibiting a total of more than 

2 x 10 7 MC sweeps for each value of q, and combining 
the results from different temperatures using the multi- 
histogram approach 0, see also Ref. 0. For q — 4.05 a 
two-peak structure is visible, while for q = 4.00 not. This 
confirms within the given numerical accuracy the known 
result q c = 4. 

For three dimensions, the situation is less clear, no ex- 
act analytic results are available. A value of 2 < g c (3) < 

3 seems likely, see analytical wo rk ll^llHll7l and simula- 
tions of moderate-size systems |l8l Il9l l20j . In the range 



FIG. 3: Free energy F per spin of the three-dimensional Potts 
model (2 < q < 3) as a function of the temperature T for 
system size L = 100. The arrow indicates the phase transition 
for q = 3. The inset shows the average number of edges Nt 
as a function of T for q — 2.1, 2.4 and 2.7. 

where the transition is first order, the transition seems 
to be weak, i.e. making a direct numerical treatment 
difficult. This is confirmed by the results for the free en- 
ergy calculated using the present approach for N = 100 3 , 
see Fig. The data is obtained by combining the re- 
sults for k — 212 different values qi G [1,3]. No clear 
kink in any of the functions is visible. Thus, to ob- 
tain a precise estimate for q c , one has to study again 
the number of edges. The average, shown in the inset 
of Fig. |3 allows to see the transition point well, but to 
infer the order of the transition is still difficult because 
of finite-size rounding of the curves. The full distribu- 
tions close to the phase transition are presented in Fig. 
0] The distributions are obtained again by performing 
simulations for several temperatures close to T C (L, q), for 
each value of q more than 2 x 10 6 MC sweeps, and com- 
bining the data using the multi-histogram approach 0. 
For q = 2.6, 2.5 one can see a clear double-peak struc- 
ture, while for q — 2.4 the distribution has only a faint 
double-peak structure. Clearly just one peak is present 
for q = 2.3. Since the depth of the minimum between the 
two peaks grows with system size |28j . it is still possible 
that q c is even below q — 2.3, but from the shape of the 
distribution at q — 2.3 this seems unlikely. This allows 
to conclude q c (d = 3) = 2.35(5) from the present results. 

This result is smaller than the result q c = 2.620(5) ob- 
tained by Gliozzi. The deviation is probably due to the 
fact that in that work much smaller system sizes N = 14 3 
where used, which shifts the value of q c up, as discussed 
above. The result q c — 2.45(10) obtained by Lee and 
Kosterlitz 0] is compatible with our result, although is 
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FIG. 4: Distributions of the number of edges at temperature 
T close to the critical value T c (L,q) for q = 2.3,2.4,2.5 and 
2.6. 



even less reliable since small sizes were used and data 
obtained at q = 3 was extrapolated to values q £ [2.7, 3]. 
Barkema and de Boer studied [lO] a model with inte- 
ger q, but mimicking the behavior of any q > 0, and 
got q c — 2.21. The results of analytical studies are scat- 
tered around the result obtained here: Kogut and Sinclair 
found 116 1 q r = 2.55 in a 1/q expansion, Nienhuis et al. 
obtained [15j q c ~ 2.2 using a real-space renormalization 
approach, while Grollau et al. got [ljj 9c ~ 2.15 within 
a Ornstein-Zernicke Approximation. 

To summarize, a new approach to calculate numeri- 
cally the partition function of g-state Potts models for 
arbitrary values q > is presented. Using a combination 
with a fast cluster algorithm large system sizes can be 
treated. The method is evaluated by performing a com- 
parison with exact analytic results for two-dimensional 
Ising models of size N — 1000 2 , a very good agreement 
is found. For the d = 2 Potts model of the same size, the 
analytically obtined critical value q c = 4 is confirmed. 
For the three-dimensional Potts model, due to the weak- 
ness of the first order transition, it is hard to infer q c 
from the data (size N = 100 3 ) for the free energy. From 
the analysis of the distribution of the number of edges in 
the generated subgraphs, q c — 2.35(5) is concluded. 

The approach to obtain the free energy should be ex- 
tensible beyond the standard g-state ferromagnetic Potts 
model, e.g. for random or diluted ferromagnets, other 
lattice types and/or higher dimensions. The method 
should work in principle also for frustrated systems, but 
here the efficient generation of the subgraphs remains an 
open problem. 
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